Attitude Control of Ornithopter Wing by Using a MIMO Active Disturbance Rejection Strategy

This work proposes a mathematical solution for the attitude control problem of an ornithopter wing. An ornithopter is an artificial bird or insect-like aerial vehicle whose flight and lift movements are produced and maintained by flapping wings. The aerodynamical drag forces responsible for the flying movements are generated by the wing attitude and torques applied to its joints. This mechanical system represents a challenging problem because its dynamics consist of MIMO nonlinear equations with couplings in the input variables. For dealing with such a mathematical model, an Active Disturbance Rejection Control-based (ADRC) method is considered. The cited control technique has been studied for almost two decades and its main characteristics are the use of an extended state observer to estimate the nonmeasurable signals of the plant and a state-feedback control law in standard form fed by that observer. However, even today, the application of the basic methodology requires the exact knowledge of the plant’s control gain which is difficult to measure in the case of systems with uncertain parameters. In addition, most of the related works apply the ADRC strategy to Single Input Single Output (SISO) plants. For MIMO systems, the control gain is represented by a square matrix of general entries but most of the reported works consider the simplified case of uncoupled inputs, in which a diagonal matrix is assumed. In this paper, an extension of the ADRC SISO strategy for MIMO systems is proposed. By adopting such a control methodology, the resulting closed-loop scheme exhibits some key advantages: (i) it is robust to parametric uncertainties; (ii) it can compensate for external disturbances and unmodeled dynamics; (iii) even for nonlinear plants, mathematical analysis using Laplace’s approach can be always used; and (iv) it can deal with system’s coupled input variables. A complete mathematical model for the dynamics of the ornithopter wing system is presented. The efficiency of the proposed control is analyzed mathematically, discussed, and illustrated via simulation results of its application in the attitude control of ornithopter wings.


Introduction
The control of flapping-wing flight vehicles has been receiving increasing attention from the research community over the last few decades [1][2][3][4][5][6][7][8][9]. In this growing field, the literature has reported remarkable progress toward the design, construction, and control of flapping-wing aerial vehicles. Since it is inspired by biological systems such as the flights of birds and insects, it has captured the attention and interest of engineers, roboticists, hobbyists, and even biologists. This interest is also due to the ability of these systems to hover and maneuver quickly [5]. In [10], the authors describe the design and fabrication of a soft robotic dragonfly whose flapping wings are powered by dense dielectric elastomer actuators (DEAs). Some of the advantages raised by the authors are that these aerial robots operated by DEAs demonstrated insect-like resilience and agility. However, a notable drawback is that these flexible actuators require upwards of 500 volts to operate. Notably, this condition presents significant practical challenges for the use of electronic controllers, as the latter require compact and efficient high-voltage power circuits. Nevertheless, many of its attractive motion properties also come at a price. For example, the airflow generated by the wings' flapping increases the aerodynamics' complexity around the flier's control surfaces of the flier [6]. Undoubtedly, these challenges/difficulties need to be taken into account in the design of the controller algorithm.
In [9], the authors were concerned with describing the mathematical model of the ornithopter using the Rayleigh-Ritz method which, according to the reported study, is more robust to wing shape changes and seems to be more suitable for optimization of profile design. Although the authors have not addressed the control problem of the mechanism, the proposed model seems to be computationally efficient and suitable for control purposes, and complex parametric analyses. In [11], the authors present an analysis of the unstable aerodynamics of low-aspect-ratio ellipsoidal wing-flapping ornithopters. In the study reported by these authors, the aerodynamics are analyzed and modeled using three-dimensional simulations by the Computational Fluid Dynamics (CFD) method.
One of the most significant difficulties in the control systems design is obtaining the mathematical description of the system since a modeling error may generate unmodeled dynamics and/or uncertainties in the values of the system's parameters. The system is then said to be uncertain, as its control is a relevant problem because of its industrial and academic applications. Among the several control strategies for uncertain systems proposed in the literature concerning linear and nonlinear systems, it can be cited the adaptive Backstepping [12], Model Reference Adaptive Control (MRAC) [13], VS-MRAC [14], and others. In addition to the system uncertainties, another fundamental problem in the control system design is the efficient rejection of external disturbance, where some relevant nonlinear control strategies were proposed [15][16][17][18].
In this context, the Active Disturbance Rejection Control (ADRC) strategy emerged [19][20][21][22], which uses an Extended State Observer (ESO) to estimate the unavailable signals of systems, including the nonmodeled dynamics and the external disturbance. Then, these estimates are used in a stabilizing control law, generating a control strategy with simple implementation and robustness properties to external disturbance and nonmodeled dynamics. Due to these advantages, extensions of this control strategy have been developed, allowing its application in a broader class of systems. As an example, in Zheng et al. [23], the ADRC strategy was extended to Multiple Inputs-Multiple Outputs (MIMO) systems being applied in an industrial chemical process. For higher-order systems, in [24], it is proposed the use of a cascaded of the reduced-order estimators, where the estimation and control parameters are maintained for each observer, and the estimates of the parametric values are obtained using a particle swarm optimization procedure according to an established cost function. In Sun et al. [25], an extension of ADRC is proposed for nonminimum phase systems, where a feedforward control assisted by the model is used to assure a null stationary error with a minimum settling time. As can be seen in these and other recent works involving ADRC for SISO plants [21,22,24,[26][27][28][29], and also for MIMO plants [30,31], the uncertainties in dynamical parameters are taking into account, despite assuming that the control gains are considered to be known. Therefore, in many practical applications, such as load transportation on mobile robots and/or drones, the values of load mass and the system gravity center vary during the task execution, which causes an uncertainty in the values of the inertia matrix. This matrix strongly influences the control gain [32,33], which also becomes uncertain.
Aiming to solve the problem of an uncertain control gain, an extension of ADRC, namely modified ADRC, was proposed in [34]. The modification consists of inserting a constant gain in series with the system's output and a filter inserted in parallel with them. Then, aiming to obtain an estimate considering the external disturbance, nonmodeled dynamics, and parametric uncertainty (including the control gain), an ESO is designed, where it is only necessary to know the signal of the control gain.
However, as only the case of SISO (Single Input, Single Output) systems was mathematically analyzed, the control gain was considered a scalar constant. Thus, to MIMO systems, the mathematical analysis of [34] could be applied only if the input was decoupled (or a weak coupling), where the systems could be represented as a set of SISO systems. However, many MIMO systems have a strong coupling between the inputs. In this case, the strategy of considering the MIMO system as a set of SISO systems has not proved to be efficient. ADRC MIMO strategies were also recently proposed in [35,36]. In [35], the authors describe a very interesting application involving trajectory tracking control in a bus-type land vehicle, in an approach called Self-Driving Bus (SDB). In the case of this system, there is a tight coupling between the vehicle lateral error subsystem and the vehicle speed subsystem. The authors then develop ADRC-based decoupling control for a class of general uncertainty MIMO nonlinear systems. In [36], the application of the ADRC method in the control of Continuous Robots (CRs) is discussed. These are a class of robotic systems with a series of connecting units gathered to perform a given task in a continuous motion manner. In this work, the CR of interest is modeled as a Multi-Input Multi-Output (MIMO) system subjected to an external disturbance, and whose internal parameters are considered uncertain. The ADRC MIMO strategy used in the latter work is an extension of the standard ADRC strategy for SISO systems in which the uncertain part of the control matrix gain is included in the generalized disturbance.
This work presents and discusses an extension of the SISO Modified-Plant ADRC (MP-ADRC) to MIMO systems. As part of the analysis, it is shown that all the robustness properties of the standard SISO MP-ADRC are kept unchanged when applied to the MIMO case. Also, an additional advantage is highlighted in what concerns the knowledge of the system's matrix control gain. Such a priori information is not required since only the signs of its main diagonal entries are needed, which contracts with other ADRC MIMO strategies. In order to illustrate the application of the proposed MIMO MP-ADRC strategy as a case study, the work considers the attitude control problem of an ornithopter wing, which is also addressed by other works in the literature. For example, the Proportional-Integral-Derivative (PID) controller and its variants are the commonly used schemes for flight control in real-world applications [7,8,37,38]. Classical optimal techniques such as Linear Quadratic Regulator (LQR), Linear Quadratic Integrator (LQI), and Model Predictive Control (MPC) are also adopted as control solutions [39,40]. Robust control schemes such as H ∞ also appear in literature as a choice. Sliding Mode Control (SMC) is another robust scheme reported in the literature [41]. Active disturbance rejection techniques are also reported in the works [31,[42][43][44].
In [3], the authors investigated the hovering performance of a Micromechanical Flying Insect [1] by designing LQR-based feedback control laws. In particular, this work provides a methodology to approximate the time-varying dynamics governed by the aerodynamic forces with a time-invariant mathematical model obtained via averaging theory and a biomimetic parametrization of the wing trajectories. In [45], the proposed controller is based on adaptive backstepping and requires no knowledge of the aerodynamical parameters. The control algorithm can follow a given angular reference trajectory during the flight task only by actuating on the tail deflection. The work reported in [46] discusses the flight stabilization problem of an avian-scale flapping robot with articulated wings. Because of the impact of the fuselage pitch angle influence on the flight altitude, a cascade control scheme is used to control fuselage and tail pitch angles in inner loops and the altitude in the outer one. In both loops, classical PID control laws are used. In the works [7,37], the authors have used feedback linearization plus PD control methods, which are the most common control approaches applied for nonlinear robotics systems.
In [42], authors use two isolated ADRC schemes to stabilize the attitude of the pitch and roll angles. In the proposed control strategy, disturbances and the coupling effects between pitch and roll channels are estimated by an Extended State Observer (ESO) and then compensated in real time. In [44], an interesting control problem is addressed. The work discusses the stabilization problem applied to flapping-wing robots before a take-off phase. The authors propose a control scheme composed of a feedback linearization part and an active Disturbance Rejection Control (ADRC) estimation stage for rejection purposes.
Although the cited ADRC-based schemes have been well documented in the literature, few have explored the control problem from the multivariable systems point of view. In this sense, the present work intends to present a valuable contribution to the state of the art. In summary, this paper intends to present the following contributions: 1.
Present a detailed extension of the MP-ADRC method for application in MIMO systems; 2.
Discuss and analyze such extension in the attitude control of an ornithopter wing, which corresponds to a relevant and challenging nonlinear control problem; 3.
Analyze the results of computational simulation obtained after control implementation, taking into account the detailed dynamic model of the ornithopter wing; 4.
Compare the performance of the proposed control scheme with those obtained with the traditional feedforward PD strategy [8,47].
This paper is organized as follows: in Section 2, the problem statement is presented. The computed torque, a traditional control strategy, is presented in Section 3. The standard ADRC and the MP-ADRC for MIMO systems, proposed in this paper, are presented in Section 4 and Section 5, respectively. Section 6 approaches the application of MP-ADRC for MIMO systems in the dynamic model of the ornithopter wing. The simulation results are shown in Section 7. Finally, Section 8 approaches the final conclusions.

Problem Statement
In this paper, we consider the attitude control problem of one ornithopter wing, as illustrated in the images of Figure 1. The mechanism of the wing model considered in the current work consists of a three-degrees-of-freedom rotary spherical joint that orientates a flat surface wing according to the joint variables α, β, γ ∈ IR [rad]. Such rotation angles are defined with respect to their corresponding axis, Z, Y, and X, respectively.

Bird Wing Insect Wing
Wing Model The wing motion is performed through the application of the input torques τ x , τ y , τ z ∈ IR [Nm] around the Z 0 , Y 0 , and X 0 axis, respectively. This way, the control objective is to derive convenient mathematical laws for the torque inputs in order to force the joint variables to track some desired angular profiles represented here by α * , β * , γ * ∈ IR [rad]. The control system general overview is shown in Figure 2, in which the output error variables e α , e β , e γ ∈ IR [rad] are highlighted together with an illustrative image of the desired angular profiles. Related to the control problem addressed here, the next section presents and discusses the Computed Torque method, which is a traditional control strategy [7,8,37,38].

Computed Torque Method
As reported in several works in the literature, e.g., [8], the ornithopter wing motion can be represented by the well-known dynamical equation: in which q = [α β γ] T is the joint vector, τ = τ α τ β τ γ T is the torque (input) vector, T is the vector of aerodynamical forces that act in the wing, M(q, q) ∈ IR 3×3 is the symmetric positive definite inertia matrix, C q (q,q) ∈ IR 3×3 is the Coriolis matrix, and N(q) ∈ IR 3 is the vector of gravity force terms and other generalized forces that act in the wing joints. All these dynamical parameters are detailed in Appendix A. Since M(q,q) is nonsingular due to its well-known structural property [47], then Equation (1) can be rewritten asq or, in the following equivalent compact form From Equation (2) to Equation (3), the description of the dynamics has been conveniently rearranged to represent the set of original equations in an equivalent (still vector) set. This equivalent one now consists of a linear double integrator plant subjected to the vector input terms F(q,q, τ v ) and B(q,q)τ, as depicted in Figure 3. In the ideal case, considering that all model parameters and variables are known and available, a control law for the torque vector τ could be chosen as follows.
in which K 1 , K 2 ∈ IR 3×3 are diagonal and positive definite matrices with convenient entries. By simply replacing the expression of Equation (4) in Equation (3), one can verify that the closed-loop system assumes the format: which represents an exponentially stable MIMO error system. However, in many applications, the entries of the vector term F(q,q, τ v ) can involve uncertain parameters, unknown and/or unmodeled dynamics, unmeasurable signals, and external disturbances, which prevent Equation (4) from being implemented directly as it is. For solving the MIMO control problem concerning ornithopter dynamics, we adopt an extension of the SISO ADRC method proposed in [34]. The main characteristics that make ADRC an attractive control technique for the current challenge are its well-known robustness properties against parametric uncertainties, the unmodeled dynamics of the system's model, and its ability to reject external disturbances. In the basic SISO ADRC method introduced in [19], and further improved in [21,48,49], the design procedures consider the knowledge of the system's control coefficient (or simply the control gain) or some nominal value nearby. In [34], such a priori knowledge is avoided, and only the sign of that coefficient is required in the control design. Although the latter method has been proposed for general SISO nonlinear systems, this work exploits an extension for MIMO systems aiming to take advantage of some of its established robustness properties. An illustrative diagram of the proposed MIMO ADRC scheme is also shown in Figure 3.

The ADRC Method Applied to MIMO Dynamical Systems
Let us consider a class of nonlinear MIMO plants in the following format.
in which X = [x 1 x 2 x 3 ] T ∈ IR 3 represents the vector of the system's output variables, 3 represents the vector of the system's input (control) variables, 3 represents the vector of the system's nonlinear functions, and W = [w 1 w 2 w 3 ] T ∈ IR 3 is the vector of external disturbances. For illustrative purposes, we can expand Equation (8) to highlight its line expressions: in which b ij (i, j = 1, 2, 3) are the constant coefficients for the control variables τ 1 , The control objective is to define convenient laws for τ 1 , τ 2 , τ 3 that force the system's output variables x 1 , x 2 , x 3 to converge asymptotically to the desired bounded reference trajectories x * 1 , x * 2 , x * 3 , respectively. For simplifying notations in Equation (8), the vector F(X,Ẋ, W) will be denoted by F(t) henceforth.
Regarding the general MIMO plant of Equation (8), from Equation (4), it is well-known that for a nonsingular B and known F(t), a control law that would achieve the desired output tracking ideally would be given by: where T is the vector of reference trajectories and K 1 , K 2 ∈ IR 3×3 are diagonal and positive definite matrix gains. However, in many applications, the generalized disturbance vector F(t) can involve unknown and/or unmodeled dynamics, unmeasurable signals, and external disturbances. In such cases, the control law of Equation (10) can not be directly computed.

Standard ADRC Design Methodology
In order to overcome the implementation problem of Equation (10), an Extended State Observer (ESO) is designed to estimate F(t) and all the system state variables. The strategy is to define the plant state vector asX and insert F(t) as an additional state variable: By assuming that F(t) is a vector of differentiable functions, the system state-space representation assumes the following format:X in which 0 ∈ IR 3×3 is the zero matrix and I ∈ IR 3×3 is the identity matrix. Note that in such a configuration, the pair (A, C) is always observable. Thus, the design of an Extended State Observer (ESO), with stateX is feasible and can be given by: T is the estimate ofX and L ∈ IR 9×3 is the estimator gain matrix. Note, from Equations (12) and (13), that estimation error dynamicX =X −X is described by˙X Moreover, since Equation (12) is an observable representation, the eigenvalues of A − LC can be placed in predetermined locations on the Left Half-Plane (LHP). In this situation, the system represented by Equation (14) is BIBO, and the estimation errorX will remain bounded for a boundedḞ(t), ∀ t. Note also that due to the property of pole placement freedom, the error convergenceX → 0 will be guaranteed in the caseḞ(t) tends to zero. Thus, using the estimated statesX in Equation (10), the following state-feedback control law is proposed: Stability Analysis Let us define the reference trajectory vector asX * := [X * T ,Ẋ * T , 0] T and the tracking error asĒ :=X −X * . Then, from Equation (12), the tracking error dynamics can be described by˙Ē Note that the control law τ can be defined like Equation (15), i.e., where we define K := [K 1 , K 2 , I]. Then, by replacing Equation (17) into Equation (16), from which we highlight the following propertȳ Thus, by gathering together Equation (14) and the remaining part of Equation (18) as an augmented error vector, then the overall system closed-loop error dynamics results in: Since the closed-loop characteristic roots emerge from the eigenvalues of matrices A − LC and A −BK, these can be placed arbitrarily on the LHP by choosing the entries of both the ESO gain matrix L and the controller gain matrix K conveniently. Hence, by assuming a boundedḞ, which is a reasonable assumption in practice, the following statements hold: • The resultant closed-loop system is BIBO; • The tracking error entries of vectorĒ are uniformly bounded ∀ t; • In a steady-state regime,Ē → 0 ifḞ(t) → 0.
In the previous approach, note that B is required to be a known matrix for implementing the control strategy of Equation (15). If B is not known exactly, then it is not possible to implement the ESO of Equation (12) as it is. Thus, the generalized disturbance F(t) estimate can not be guaranteed to succeed, leading to tracking error degradation or instability. The problem concerning the lack of requirement on the knowledge of the B matrix is addressed in the next section, in which an extension of the MP-ADRC method is discussed for MIMO systems with uncertain B matrix.

Modified-Plant ADRC (MP-ADRC)
Consider the general class of dynamical systems described by Equation (8). For analysis purposes, let us detach the system linear apart from F(t) so that it can appear explicitly in the plant representation, as follows: where A 1 , A 0 ∈ IR 3×3 are constant matrices and G(t) ∈ IR 3 is the generalized disturbance vector that includes the nonlinear dynamics and the external disturbances. Before continuing to develop the MP-ADRC, the following assumptions are considered: Assumption 1. The entries g i (t), (i = 1, 2, 3) of the generalized disturbance vector G(t) and their first-order derivativesġ i (t) are bounded by known constants g p , g d > 0 ∈ IR, that is, Assumption 2. The constant entries of A 1 , A 0 , and B are uncertain, but with known constant upper and lower bounds M 0 > 0 ∈ IR and M b > 0 ∈ IR, respectively, given by: where σ m (.), and σ M (.) stand for the lowest and the largest matrix singular value, respectively.

Assumption 3.
The matrix B is nonsingular, dominant diagonal and the signs of its main diagonal entries are known.
As will be shown in the next section, the Assumption 1 is critical for choosing the bandwidth of the state estimator, which should be greater than the bandwidth of G(t). Assumption 2 allows us to consider a more general class of plants with a complete set of uncertain parameters, including the input coefficient (also known as control gain). Assumption 3 is feasible since several mechanical systems have at least a main actuation on each degree of freedom.

Proposed Methodology
The idea of the modified ADRC scheme is to perform a structural change in the plant of Equation (8) input/output description in order to obtain a new dynamical system with some desired characteristics. As will be shown, by adopting such a modification, both controller and ESO designs are simplified in terms of fewer requirements on the knowledge of plant parameters. The original system tracking objective is kept unchanged together with the well-known stability and convergence properties.
As can be seen in the diagram of Figure 4, the plant modification is performed by a constant gain matrix β o ∈ IR 3×3 that is introduced in series with the plant output error and a diagonal matrix of second-order filters D(s) in parallel with them. These filters are given by: in which, is a stable filter (poles with negative real parts). Based on the diagram of Figure 4, it can be concluded that: where Λ 1 = α 1 I and Λ 0 = α 0 I. Note from Equation (26) that α 1 = 2p and α 0 = p 2 . By differentiating Equation (27) twice with respect to time and using the expressions in Equations (21) and (28), one can conclude thaẗ Since τ f = Z − β o E, then, from Equation (29), we have thaẗ Thus, by replacing the expression of Equation (31) in Equation (30), we obtain the modified plant equation, now written in terms of the new output error vector Z, the new generalized disturbance vector Ω(t) and the new control vectorτ: Note, from Equations (32) and (33), that the input matrix B is inserted into the new generalized disturbance Ω(t) and the new input signal isτ, which has the identity matrix I as its coefficient. Therefore, if vector Ω(t) was known, then it would be enough to choose the trajectory tracking control law asτ = −Ω(t). This particular choice could make the closed-loop system of Equation (32) become exponentially stable. However, as Ω(t) is an uncertain vector, then an ESO is suggested for estimating it. Thus, definingZ 1 = Z,Z 2 =Ż andZ 3 = Ω(t) as state variables, the Equation (32) can be describe by  (34) does not require the value of the input gain matrix B of the plant, which contrasts with other ADRC methodologies in the literature. This is due to the fact that now B is inserted within the disturbance Ω, as described by Equation (33). In this approach, B will be estimated together with Ω. This is the main feature that gives the MP-ADRC method a higher level of robustness to uncertainties in B.

Remark 1. Note that ESO's implementation of Equation
In order to facilitate the estimator design, a state transformation matrix T is used: is the state vector in the canonical observable form. Note that Z =Z 1 =Z o1 and Ω(t) =Z 3 =Z o3 . Then, it is known thaṫZ whereĀ is the ESO gain matrix that is responsible for allocating the closed-loop characteristics roots (or poles) of the estimator. The entries ofL can be computed by solving the matrix polynomial equation being ω 0 ∈ IR a positive design constant that defines the ESO characteristics roots and also defines the size of its bandwidth. Then, one can conclude that Hence, the control law to the modified plant of Equation (32) can be described now as a function ofẐ o3 , that is the estimate of Ω(t): Note that the convergence analysis ofẐ o3 to generalized disturbance Ω(t) is essential to verify the stability of the closed loop. This analysis is performed in the next section.

ESO Convergence Analysis
From Equations (36) and (37), it can be concluded that the estimation error equation is described by˙Z From Equation (46), one can conclude that ...
Therefore, in the frequency domain, from Equations (40)-(42), we have that As can be seen from Equation (44), E Ω = Ω −Ẑ o3 and then, after a few algebraic manipulations, the expression in Equation (47) in which emerges the scalar transfer function H(s) that represents the dynamical relation between the entries of the generalized disturbance Ω(t) and their estimates withinẐ o3 . In the frequency domain, due to the low-pass behavior of H(s) in Equation (48), its modulus always lies on the interval (0; 1]. This means that if ω 0 is chosen with a sufficiently large value greater than the highest frequency of Ω(t), then Equation (48) can be approximated byẐ o3 (t) ≈ Ω(t). (47) and (48) is only for analysis purposes. Moreover, in cases where Ω(t) involves nonlinear functions of the plant state variables, the definition of an Ω(s) in the frequency domain may not be consistent. Thus, as we want to analyze the amplitude of the estimation error signal in Equation (47) without losing generality, we believe that, in this work, this is a more appropriate mathematical formalism to represent the input/output relationships.

Closed-Loop Stability and Tracking Analysis
By assuming that ω 0 in Equation (48) is chosen conveniently, let us define the variable c 0 ∈ (0; 1] for representing the value of the equivalent gain of |H(jω)| computed in some frequency ω, [rad/s] within the range [0; ω 0 ). Such a definition is introduced here only for analysis since we aim to focus on the magnitude behavior of the error system in a steady-state regime. For this purpose, we have that c 0 = |H(jω)|. (49) Thus, based on Equation (49), we can rewrite the control law in Equation (43) aṡ By replacing Equation (33) in Equation (50), we have thaṫ Therefore, using Equations (21) and (29) it can be concluded thaṫ From Equation (26), the control law (52) can be described, in frequency domain, by Note, from Equations (53) and (21), that the closed-loop system can be represented by Figure 5, where Λ = pI. By rearranging the blocks in Figure 5 to place the signal X(t) as output and the signals G(t) and X * (t) as inputs, the resulting block diagram assumes the format of that one depicted in Figure 6, where In this case, the closed-loop system will be asymptotically stable. Moreover, it can be concluded, from Equation (57) and Figure 3, that Then, by replacing Equation (57) in Equation (58), we have that Therefore, Thus, one can conclude, considering the Assumption 1, that ||E(t)|| tends to a residual set which can be reduced with the increase of ||c 0 Bβ o || ∞ . In addition, E(t) → 0 if G(t), X * (t) tend to be constant values. The following theorem summarizes the obtained result.

Remark 3.
The mathematical developments that result in Theorem 1 were applied in a MIMO system composed of three second-order equations. This system was chosen to facilitate the analysis understanding. However, all the steps can be easily extended to systems with a higher number of equations. It is well known that many physical MIMO systems can be mathematically described by n second-order equations, which shows the generality of the results discussed in the present work. For MIMO systems of general order, the design procedures of the modified ADRC technique can be suitably applied; however, they will require some adjustments in the stability-proof mechanism. This will be left as a proposal for future work.

Remark 4.
The occurrence of delays in the closed-loop system is a practical reality in control systems. As a general consensus in the literature, delay is considered an unmodeled dynamic. Although we have not emphasized this in the current analysis, from a theoretical point of view, one of the characteristics of the ADRC method is the ability to deal with this type of dynamical characteristic. We believe, to the best of our knowledge, that the developments around this analysis can be driven by using the theory presented here in this work together with the Padé approximation Function Theory. In fact, this issue needs to be investigated more closely, which could lead to extensive mathematical developments, additional simulations, and experiments. This topic is not addressed in the current work, by believing that the presentation of such an analysis is not trivial.
In this way, we judge it will be more convenient to deal with such a rigorous and complete analysis in future work.

Application-Attitude Control of an Ornithopter Wing
In this section, the modified ADRC technique, discussed in the previous sections, is applied to the attitude control problem of an ornithopter wing with a dynamical model described by Equation (1).

System Model
Since the mechanism's motions behave like a spherical joint, its attitude control can be defined as a tracking problem for the output angles α, β, γ. For the wing motion description, three parameters are introduced: the mass m, the chord C a , and the span L a , as depicted in Figure 7. The dynamical model for the wing motion is described by (1) (see Appendix A for more details).

Control Definition
Since M(q,q) is a nonsingular matrix, then Equation (1) can be rewritten as in which we definē and then,  By noticing that the ornithopter wing dynamics of Equation (63) is similar to Equation (9), with B = M −1 , then the proposed MP-ADRC scheme of the previous sections can be applied here directly. In fact, after choosing the design parameters p of Equation (26), β o (28),L (40)- (42) properly, one can directly implement Equation (29) for filters, (37) for the ESO and (43) for the control law τ. In this case, the tracking error is defined as in which the reference trajectory is now given in terms of the desired angle vector q d , namely Thus, based on the developments of the previous sections and on the similarity of systems in Equations (63) and (9), we highlight in the following the final expressions, in expanded form, for summarizing the proposed control system.
For the filter: For the auxiliary error: For the ESOs: (a) Subsystem α:˙Ẑ (c) Subsystem γ:˙Ẑ For the control laws: In the next section, numerical simulations are performed to illustrate the application of the proposed strategy. The obtained results are presented and discussed after comparisons with other classical control methods.

Simulation Results
In this section, the performance of the modified ADRC of Equations (64)-(75) is compared with a nonlinear version of the augmented PD control strategy [47] described by In addition to the equations explained in the previous section, Figure 8 shows a compact-form block diagram that illustrates the implementation of the proposed method. The comparison is performed by changing some dynamical parameter values of the ornithopter wing model to mimic some sorts of parametric uncertainties. This is performed in two different simulation trials, as shown in Table 1. The complete set of model parameters and matrices can be found in Appendix A. For comparison purposes, the values of the control parameters for both control strategies are held the same in the two simulations, as shown in Table 2. This is for allowing the robustness analysis of the control strategies. It is important to stress that these control parameter values are applied in the three-degrees-of-freedom model of the ornithopter wing. The design parameters of Table 2 were computed based on the tuning algorithms presented in references [34,51]. As argued in [51], such parameters play an important role in the stability of the closed-loop system and then need to be tuned properly. It has been also shown in [51] that the tuning procedure can follow the well-known root locus method. In both simulations, a wind velocity is considered as an input disturbance (see Figure 9), and white noise is introduced as a measurement noise ( Figure 10). Since it is known that the output derivative calculation is not an appropriate measurement for practical applications (because of inevitable measurement noise), the following first-order filter is used for computingq andq d : The use of a filtered derivative is imperative for the practical implementation of the PD algorithm since an exact derivative is not implementable in practical applications. However, this is not required in the current ADRC MIMO proposal, since the system states are estimated by the ESO. In addition, actuator saturation values of ±30 N·m are considered.  Figures 11 and 12 show the three wing joint angles (α, β, γ) and their reference trajectories (α r , β r , γ r ) considering the modified ADRC and augmented PD with a time constant τ d = 10 −4 . The tracking errors are shown for both controllers in Figures 13 and 14. Note that, by analyzing this set of four figures, the modified ADRC and augmented PD showed similar performances over time. However, as shown in Figures 15 and 16, the augmented PD presents a very noisy control signal (axis torques) when compared with the modified ADRC approach. Such behavior is due to the derivative term acting as an amplifier for the measurement noise, which is a known drawback of augmented PD. Note that the tracking error amplitude is greater than those shown in Figures 13 and 14. Moreover, although the control signal has a minor noise amplitude, its performance is worst in the initial instants, also presenting a switching behavior due to the saturation levels of ±30 N.m. As shown in Figures 19 and 20, the augmented PD performance is degraded even more when τ d = 10 −2 , which is expected due to the worse approximation of the derivative term.

Simulation 2
Figures 21-26 show the simulation results when some wing parameters are changed (see Table 1). Note that both strategies showed good robustness to the variation of dynamical parameters. However, as in simulation 1, the modified ADRC again presented a smaller tracking error amplitude and a resulting control signal with less noise.
Therefore, the simulation results show that the modified ADRC presented a similar tracking error performance as the augmented PD with τ d = 10 −4 , but the augmented PD control signal has minor robustness to measurement noise, presenting a significantly noisier control signal. By using τ d = 10 −3 , the augmented PD control signal noise is reduced, although it was also higher than that of the modified ADRC. However, the tracking error performance becomes more significant in terms of small amplitude.
Thus, besides its simplicity, since the augmented PD strategy needs a complex mathematical description of the controlled system, the modified ADRC presented a better performance, mainly due to its better robustness to measurement noise, presenting a less noisy control signal.

Conclusions
This paper discusses a design procedure for the control of MIMO uncertain systems as an extension of the modified ADRC strategy for SISO plants. The proposed strategy, which, in the SISO case, is known for its robustness properties against parametric uncertainties, unmodeled dynamics, and external disturbances, also revealed good performance for MIMO plants with coupled input variables. Developments show that the only information required is the sign of the main diagonal entries of the control matrix. The relevance of the MP-ADRC method compared to LADRC lies in the fact that the former does not depend on the control gain value B of the plant. Even in the LADRC strategy, decoupled control gains are required both for the ESO design and for the law of formation of the control variable. In the MIMO MP-ADRC, although the approach considers a matrix gain B, this is not required either in the ESO project or in the definition of the control law, as can be verified by Equations (34) and (43). In order to illustrate the efficiency of the proposed strategy, the methodology was applied to the wing attitude control problem of an ornithopter. Then, simulation trials were performed and compared with the extended PD scheme. After analyzing the obtained results, it was observed that both controllers' performance was similar regarding tracking error and robustness to uncertain parameters. However, the extended PD presented a noisier control signal due to the derivative effect. The use of filtered derivatives is imperative for the implementation of the PD algorithm. Moreover, unlike the ADRC, the robustness of the PD control proved to be dependent on τ d in the Equation (77), as can be observed from Figures 19 and 20, where the performance has deteriorated as τ d increased. On the other hand, the control signal is noisier for small τ d , as can be observed from Figure 16, which is very critical in practical applications. Therefore, the simulation results showed that the MP-ADRC has the advantage of good robustness with a less noisy control signal when compared with augmented PD control. Then, in addition to being simpler to implement, it is possible to conclude, from the simulation results, that the MP-ADRC MIMO showed a better general performance after analyzing the tracking error and the control signal.  Data Availability Statement: Data sharing is not applicable to this article.

Acknowledgments:
The authors would like to thank CEFET/RJ, the Rio de Janeiro research agency FAPERJ, and the Brazilian federal research agencies, CAPES and CNPq, for supporting this research.

Conflicts of Interest:
The authors declare no conflict of interest.